#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <limits.h>
#include <string.h>
#include <assert.h>
#include <float.h>

#include "pdist.h"

/**
 * Print a matrix to stdout.
 *
 * Primarily used for debugging purposes.
 */
void print_matrix(float *input, const unsigned int size_x, const unsigned int size_y)
{
	for (unsigned int j = 0; j < size_y; j++)
	{
		printf("%5d: ", j);
		for (unsigned int i = 0; i < size_x; i++)
			printf("%-03.3f ", input[j*size_x + i]);
		printf("\n");
	}
	printf("\n\n");
}

/**
 * Write a matrix to a file.
 *
 * Primarily used for debugging purposes.
 */
void write_matrix(const char *filename, float *input, const unsigned int size_x, const unsigned int size_y)
{
	FILE *f;
	char line[1024];

	f = fopen(filename, "w");

	for (unsigned int j = 0; j < size_y; j++)
	{
		line[0] = '\0';
		for (unsigned int i = 0; i < size_x; i++)
			sprintf(line, "%s%s%.3f", line, (i == 0 ? "" : " "), input[j*size_x+i]);
		sprintf(line, "%s\n", line);
		fputs(line, f);
	}

	fclose(f);
}

#define EPSILON 0.001

/**
 * Compare two matrixes.
 *
 * Primarily used for testing purposes.
 */
unsigned int compare_matrix(float *input, float *input2, const unsigned int size_x, const unsigned int size_y)
{
	int ret = 0;
	for (unsigned int j = 0; j < size_y; j++)
	{
		for (unsigned int i = 0; i < size_x; i++)
			if (abs(input[j*size_x + i] - input2[j*size_x + i]) > EPSILON)
			{
				printf("incorrect: %d %d %3.3f %3.3f\n", i, j, input[j*size_x + i], input2[j*size_x + i]);
				ret++;
				//return ret;
			}
	}
	return ret;
}


/**
 * Compare two matrixes using a fuzzy approach.
 *
 * Because linkage can have different answers if distances are close, we have to compare fuzzily.
 */
unsigned int
compare_matrix_linkage(float *input, float *input2, const unsigned int size_x, const unsigned int size_y)
{
	int ret = 0;
	float alist[50], blist[50];
	int alen = 0, blen = 0;

	for (unsigned int j = 0; j < size_y; j++)
	{
		for (unsigned int i = 0; i < size_x; i++)
			if (abs(input[j*size_x + i] - input2[j*size_x + i]) > EPSILON)
			{
				int found = 0;

				// look for the first number in the blist
				for (int k = 0; k < blen; k++)
				{
					if (found)
					{
						blist[k-1] = blist[k];
					}
					else
					{
						if (abs(input[j*size_x + i] - blist[k]) <= EPSILON)
						{
							found = 1;
						}
					}
				}

				// if found, we can decrease b's len
				if (found)
				{
					blen--;
				}
				// otherwise, we add the number to a's list
				else
				{
					if (alen == 49) // past length, return immediately
						return 1;
					else
						alist[alen++] = input[j*size_x+i];
				}

				found = 0;

				// look for the second number in the alist
				for (int k = 0; k < alen; k++)
				{
					if (found)
					{
						alist[k-1] = alist[k];
					}
					else
					{
						if (abs(input2[j*size_x + i] - alist[k]) <= EPSILON)
							found = 1;
					}
				}

				// if found, we can decrease a's len
				if (found)
				{
					alen--;
				}
				else
				{
					if (blen == 49)
						return 1;
					else
						blist[blen++] = input2[j*size_x+i];
				}
			}
	}

	if (alen > 0 || blen > 0)
		return 1;

	return ret;
}

